## [1] "Run Completed at 2015-12-14 21:25:30"
Which birds use which plants? What are the mechanisms that determine niche overlap in this assemblage.
We will fit several difference types of models and ask how they inform our understanding. For each model, we try to calculate a measure of discrepancy and assess model fit.
Birds use flowers that have similar length of bill shape and we expect the importance of trait matching to decrease as available resources increase.
## character(0)
## [1] "Missing Trait Information: sp"
## [2] "Missing Trait Information: Burmeistera cylindrocarpa"
## [3] "Missing Trait Information: Burmeistera multiflora"
## [4] "Missing Trait Information: Calathea ischnosiphonoides"
## [5] "Missing Trait Information: Canna jaegeriana"
## [6] "Missing Trait Information: Cavendishia tarapotana"
## [7] "Missing Trait Information: Drymonia brochidodroma"
## [8] "Missing Trait Information: Drymonia collegarum"
## [9] "Missing Trait Information: Gasteranthus pansamalanus"
## [10] "Missing Trait Information: Gasternathus lateralis"
## [11] "Missing Trait Information: Guzmania lehmanniana"
## [12] "Missing Trait Information: Guzmania squarrosa"
## [13] "Missing Trait Information: Justicia secunda"
## [14] "Missing Trait Information: Macleania recumbens"
## [15] "Missing Trait Information: Palicourea sodiroi"
## [16] "Missing Trait Information: Pitcairnia nigra"
## [17] "Missing Trait Information: Psamisia ???"
## [18] "Missing Trait Information: Psammisia ecuadorensis"
## [19] "Missing Trait Information: Stromantheus stromanthe"
## [20] "Missing Trait Information: Zingibiraceae sp."
Create a binary variable whether each observation was in a low elevation or high elevation transect. We have some species that just occur at the top of the gradient, and are not present in the sampling window of flowers at the low elevation.
Accounting for non-availability. We have to figure out which plants were sampled in which periods, and if it was sampled, the non-detection are 0 if it wasn’t the non-detection are NA. then remove all the Na’s.
## Hummingbird Low m High Index
## 1 Green-crowned Woodnymph 1359 1376.645 1390 1
## 2 Rufous-tailed Hummingbird 1368 1476.333 1380 1
## 3 Andean Emerald 1378 1386.261 1380 1
## 4 White-necked Jacobin 1368 1397.500 1422 1
## 5 Stripe-throated Hermit 1368 1440.887 1450 1
## 6 White-whiskered Hermit 1340 1412.391 1450 1
What is the temporal distribution of resources?
Appears to be a log normal distribution with the existance of rare, but intense, pulses.
## Source: local data frame [1,301 x 7]
## Groups: Bird, Plant, ID [1204]
##
## Bird Plant ID DateP Yobs Elev Transect_R
## (int) (int) (chr) (fctr) (int) (dbl) (fctr)
## 1 15 85 FH709 2014-02-27 21 1990 NA
## 2 22 93 FH616 2014-05-07 17 2450 NA
## 3 17 85 FH414 2014-04-22 13 1990 NA
## 4 24 101 FL083 2013-07-29 13 2350 NA
## 5 6 98 FL061 2013-07-10 12 1325 NA
## 6 15 18 FH303 2013-11-19 12 1940 NA
## 7 15 85 FH1213 2015-02-11 11 1990 NA
## 8 24 101 FL083 2013-07-28 11 2350 NA
## 9 5 3 FL057 2013-07-05 10 1390 NA
## 10 8 92 NF084 2014-05-24 10 1513 NA
## .. ... ... ... ... ... ... ...
What elevation transect is each observation in? The camera data need to be inferred from the GPS point.
Before we account for elevation and availability, what is the distribution of intensity of interactions?
If we tend to think of the intensity of events as poisson, our data is extremely ones-inflated. The first moment of the poisson does not well capture the data.
The best fitting moment is in opaque red. The observed data is in black.
What about negative binomial?
The best fitting moment is in opaque red. The observed data is in black.
What about log normal?
The best fitting moment is in opaque red. The observed data is in black.
We have more information than just the presences, given species elevation ranges, we have absences as well. Absences are birds that occur at the elevation of the plant sample, but were not recorded feeding on the flower.
Remerge data with predictors
Three measures of resources.
Give binary versions of these resources as well - what matters is if its high or low for that resource measure. This helps account for the enormously long tail of resource pulses and minimizes the noise associated with such a widely ranging variable. The points on the huge blooms are going to have large leverage.
View quantile thresholds (0.5,.75,.9) as vertical lines.
Do we have enough data?
Relationship among resource measures
How many data points for each species for each resource level?
Let’s look at the same figures as above, with the addition of informed 0’s for non-detections.
Still not quite understanding the difference, what plants are used for each species that are not in common.
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1] "~$tworkTime.docx" "Bayesian"
## [3] "D3.R" "D3Network.html"
## [5] "figure" "Figures"
## [7] "humbipartite.html" "humbipartiteForce.html"
## [9] "humbipartiteForceNew.html" "InputData"
## [11] "Network.R" "NetworkData.Rdata"
## [13] "NetworkSource.R" "NetworkTime.html"
## [15] "NetworkTime.RData" "NetworkTime.Rmd"
## [17] "NetworkTime.Rproj" "NetworkTime_cache"
## [19] "Observed.html" "Observed.RData"
## [21] "Observed.Rmd" "PoissonSimulation.Rmd"
## [23] "README.md"
For initial inference, let’s look at some reasonable mixed effect models. There is no zero inflation here, so we treat these data with some caution.
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch + (1 | Bird)
## Data: obsf
##
## AIC BIC logLik deviance df.resid
## 13991.3 14014.0 -6992.6 13985.3 14624
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.943 -0.441 -0.286 -0.164 101.435
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bird (Intercept) 1.064 1.031
## Number of obs: 14627, groups: Bird, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.41516 0.27812 -5.088 3.61e-07 ***
## Traitmatch -0.71907 0.02997 -23.993 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Traitmatch -0.077
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch * BAll_Flowers + (1 | Bird)
## Data: obsf
##
## AIC BIC logLik deviance df.resid
## 13987.1 14025.1 -6988.6 13977.1 14622
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.982 -0.440 -0.282 -0.163 102.716
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bird (Intercept) 1.056 1.028
## Number of obs: 14627, groups: Bird, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.33549 0.32624 -4.094 4.25e-05 ***
## Traitmatch -0.51357 0.14134 -3.633 0.00028 ***
## BAll_Flowers -0.08316 0.17615 -0.472 0.63684
## Traitmatch:BAll_Flowers -0.21354 0.14380 -1.485 0.13755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc BAll_F
## Traitmatch -0.419
## BAll_Flowrs -0.524 0.773
## Trtmtc:BA_F 0.413 -0.977 -0.787
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch * BUsed_Flowers + (1 | Bird)
## Data: obsf
##
## AIC BIC logLik deviance df.resid
## 13992 14030 -6991 13982 14622
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.984 -0.443 -0.280 -0.162 100.457
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bird (Intercept) 1.069 1.034
## Number of obs: 14627, groups: Bird, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.44873 0.27907 -5.191 2.09e-07 ***
## Traitmatch -0.69960 0.03368 -20.775 < 2e-16 ***
## BUsed_Flowers 0.11819 0.06532 1.810 0.0704 .
## Traitmatch:BUsed_Flowers -0.07101 0.06138 -1.157 0.2473
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc BUsd_F
## Traitmatch -0.091
## BUsed_Flwrs -0.068 0.345
## Trtmtc:BU_F 0.045 -0.460 -0.686
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch * BFlowerA + (1 | Bird)
## Data: obsf
##
## AIC BIC logLik deviance df.resid
## 13926.8 13964.7 -6958.4 13916.8 14622
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.018 -0.454 -0.287 -0.165 92.802
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bird (Intercept) 1.124 1.06
## Number of obs: 14627, groups: Bird, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.38216 0.28662 -4.822 1.42e-06 ***
## Traitmatch -0.66934 0.03122 -21.439 < 2e-16 ***
## BFlowerA -0.08594 0.06780 -1.267 0.205
## Traitmatch:BFlowerA -0.35733 0.07352 -4.861 1.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc BFlwrA
## Traitmatch -0.085
## BFlowerA -0.058 0.312
## Trtmtch:BFA 0.034 -0.345 -0.664
## Data: obsf
## Models:
## tr_interactionA: Yobs ~ Traitmatch * BAll_Flowers + (1 | Bird)
## tr_interactionU: Yobs ~ Traitmatch * BUsed_Flowers + (1 | Bird)
## tr_interactionF: Yobs ~ Traitmatch * BFlowerA + (1 | Bird)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## tr_interactionA 5 13987 14025 -6988.6 13977
## tr_interactionU 5 13992 14030 -6991.0 13982 0.000 0 1
## tr_interactionF 5 13927 13965 -6958.4 13917 65.237 0 <2e-16
##
## tr_interactionA
## tr_interactionU
## tr_interactionF ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
What does the mixed effect model look like?
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(Yobs) ~ Traitmatch * FlowerA + (1 | Bird)
## Data: obsf[obsf$Yobs > 0, ]
##
## REML criterion at convergence: 2338.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0389 -0.6771 -0.5438 0.5315 4.6721
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bird (Intercept) 0.006509 0.08068
## Residual 0.339359 0.58255
## Number of obs: 1296, groups: Bird, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.972e-01 3.473e-02 11.439
## Traitmatch -2.177e-02 2.204e-02 -0.988
## FlowerA 1.685e-05 1.183e-05 1.424
## Traitmatch:FlowerA -2.136e-05 1.411e-05 -1.514
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc FlowrA
## Traitmatch -0.525
## FlowerA -0.259 0.249
## Trtmtch:FlA 0.126 -0.235 -0.608
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
\[ Yobs_{i,j,k} \sim Binom(Detect_i,N_{i,j,k}) \]
\[ N_{i,j,k} \sim Pois(\lambda_{i,j,k}) \]
\[log(\lambda_{i,j,k})<-\alpha_i + \beta_i * abs(Bill_i - Corolla_j) * \beta_2 * Resources_k + \beta_3 * Resources_k * abs(Bill_i - Corolla_j) \]
Priors
\[ Detect \sim N(dprior,tau_detect) \] \[\alpha_i \sim N(intercept,\tau_{\alpha})\] \[\beta_{1,i} \sim N(\gamma_1i,\tau_{\beta_1})\] \[\beta_{2,i} \sim N(\gamma_2i,\tau_{\beta_2})\] \[\beta_{3,i} \sim N(\gamma_3i,\tau_{\beta_3})\]
Hyperpriors
Group Level Means \[dprior \sim U(0,0.5)\] \[\gamma_{1,i} \sim N(0,0.0001)\] \[\gamma_{2,i} \sim N(0,0.0001)\] \[\gamma_{3,i} \sim N(0,0.0001)\] \[ intercept \sim N(0,0.0001)\]
Group Level Variance
\[\tau_{\alpha} \sim Gamma(0.0001,0.0001)\] \[\tau_detect \sim Gamma(0.0001,0.0001)\] \[\tau_\beta1 \sim Gamma(0.0001,0.0001)\] \[\tau_\beta2 \sim Gamma(0.0001,0.0001)\] \[\tau_\beta3 \sim Gamma(0.0001,0.0001)\]
Derived quantities
\[\sigma_{int} = \sqrt[2]{\frac{1}{\tau_\alpha}}\] \[\sigma_{detect} = \sqrt[2]{\frac{1}{\tau_{detect}}}\] \[\sigma_{slope1} = \sqrt[2]{\frac{1}{\tau_\beta1}}\] \[\sigma_{slope2} = \sqrt[2]{\frac{1}{\tau_\beta2}}\] \[\sigma_{slope3} = \sqrt[2]{\frac{1}{\tau_\beta3}}\]
##
## sink("Bayesian/NmixturePoissonRagged.jags")
##
## cat("
## model {
## for (x in 1:Nobs){
##
## # Covariates for true state
## log(lambda[Bird[x],Plant[x],Time[x]]) <- alpha[Bird[x]] + beta1[Bird[x]] * traitmatch[x] + beta2[Bird[x]] * resources[x] + beta3[Bird[x]] * resources[x] * traitmatch[x]
##
## #True State
## N[x] ~ dpois(lambda[Bird[x],Plant[x],Time[x]] )
##
## #Observation Process
## Yobs[x] ~ dbin(detect[Bird[x]],N[x])
##
## #Assess Model Fit
##
## #Fit discrepancy statistics
## eval[x]<-detect[Bird[x]]*N[x]
## E[x]<-pow((Yobs[x]-eval[x]),2)/(eval[x]+0.5)
##
## ynew[x]~dbin(detect[Bird[x]],N[x])
## E.new[x]<-pow((ynew[x]-eval[x]),2)/(eval[x]+0.5)
##
## }
##
## for (i in 1:Birds){
## logit(detect[i]) <- dtrans[i]
## dtrans[i] ~ dnorm(dprior,tau_detect)
## alpha[i] ~ dnorm(intercept,tau_alpha)
## beta1[i] ~ dnorm(gamma1,tau_beta1)
## beta2[i] ~ dnorm(gamma2,tau_beta2)
## beta3[i] ~ dnorm(gamma3,tau_beta3)
## }
##
##
## #Hyperpriors
## #Slope grouping
## gamma1~dnorm(0,0.0001)
## gamma2~dnorm(0,0.0001)
## gamma3~dnorm(0,0.0001)
##
## #Intercept grouping
## intercept~dnorm(0,0.0001)
##
## #detect grouping
## dprior~dnorm(0,1)
##
## # Group intercept variance
## tau_alpha ~ dgamma(0.0001,0.0001)
## sigma_int<-pow(1/tau_alpha,0.5) #Derived Quantity
##
## #Group Slope Variance
## tau_beta1 ~ dgamma(0.0001,0.0001)
## tau_beta2 ~ dgamma(0.0001,0.0001)
## tau_beta3 ~ dgamma(0.0001,0.0001)
##
## sigma_slope1<-pow(1/tau_beta1,0.5)
## sigma_slope2<-pow(1/tau_beta2,0.5)
## sigma_slope3<-pow(1/tau_beta3,0.5)
##
## #Group Detection variance
## tau_detect ~ dgamma(0.0001,0.0001)
## sigma_detect<-pow(1/tau_detect,0.5)
##
## #derived posterior check
## fit<-sum(E[]) #Discrepancy for the observed data
## fitnew<-sum(E.new[])
##
## }
## ",fill=TRUE)
##
## sink()
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph Size: 236352
##
## Initializing model
Re-update if needed to extend
Hold resource relationship at 0, what does trait-matching look like?
That looks pretty crazy, the posterior is very broad, what about holding at mean value for all coefficients
Let’s take a closer look at distribution of interaction effect posteriors values for each species.
Do species with long bill lengths have positive interaction effects?
Detection Table
## Hummingbird mean lower upper
## 5 Fawn-breasted Brilliant 0.2 0.1 0.4
## 8 Green-fronted Lancebill 0.7 0.3 1.1
## 9 Purple-bibbed Whitetip 0.7 0.3 1.9
## 10 Speckled Hummingbird 1.1 0.7 1.8
## 4 Collared Inca 1.9 1.2 2.8
## 6 Gorgeted Sunangel 2.0 1.3 3.1
## 11 Stripe-throated Hermit 2.8 1.8 4.4
## 13 Violet-tailed Sylph 2.9 2.0 4.7
## 2 Brown Inca 3.6 2.5 5.0
## 7 Green-crowned Woodnymph 4.6 1.7 12.0
## 1 Booted Racket-tail 6.5 3.5 11.3
## 12 Tawny-bellied Hermit 7.4 3.6 13.1
## 3 Buff-tailed Coronet 7.7 3.0 14.8
## 14 White-whiskered Hermit 8.1 4.1 14.2
Hierarchical Posteriors
Species level posteriors
Order predictions from worst to best.
## Hummingbird Iplant_Double chisq Yobs
## 1 Speckled Hummingbird Palicourea lineata 346.19563 17
## 2 Brown Inca Mezobromelia capituligera 222.45458 21
## 3 Gorgeted Sunangel Psammisia sodiroi 175.51353 13
## 4 Violet-tailed Sylph Mezobromelia capituligera 134.33745 13
## 5 Gorgeted Sunangel Psammisia sodiroi 131.17282 11
## 6 Brown Inca Bomarea pardina 122.64865 12
## 7 Green-fronted Lancebill Palicourea acetosoides 106.60609 8
## 8 Collared Inca Pitcairnia sodiroi 103.01062 9
## 9 Speckled Hummingbird Macleania stricta 92.34986 8
## 10 Brown Inca Heliconia impudica 86.01011 10
## Used_Flowers Traitmatch
## 1 -0.6050997 0.140
## 2 -0.5645269 0.156
## 3 1.8609052 0.408
## 4 -0.5514037 1.576
## 5 1.8609052 0.408
## 6 -0.3296283 1.988
## 7 -0.6061885 0.964
## 8 -0.5633807 1.576
## 9 0.1962717 0.536
## 10 0.3987921 1.585